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Abstract 

We study from the numerical point of view, instabilities developed in 
a fluid layer with a free surface, in a cylindrical container which is non- 
homogeneously heated from below. In particular we consider the case in 
which the applied heat is localized around the origin. An axysimmetric 
basic state appears as soon a non-zero horizontal temperature gradient 
is imposed. The basic state may bifurcate to different solutions depend- 
ing on vertical and lateral temperature gradients and on the shape of the 
heating function. We find different kinds of instabilities: extended pat- 
terns growing on the whole domain which include those known as targets 
and spirals waves. Spirals are present even for infinite Prandtl number. 
Localized structures both at the origin and at the outer part of the cyhn- 
der may appear either as Hopf or stationary bifurcations. An overview of 
the developed instabilities as a function of the dimesionless parameters is 
reported in this article. 
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Spiral and target waves instabilities in Rayleigh-Bernard convec- 
tion were experimentally found in early nineties and several theoret- 
ical mechanisms have been proposed to justify them. In this article 
we prove that these and other structures may be understood as the 
linearly growing modes of non trivial basic states generated by a non 
uniform heating function. In this framework we find that also local- 
ized structures may appear. These patterns are localized at different 
positions in the domain depending on the basic flow from which they 
bifurcate. This article reports a sytematical numerical study on these 
transitions. 
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1 Introduction 



Instabilities and pattern formation in buoyant flows have been extensively stud- 
ied in recent years. Classically, heat is applied uniformly from below and 
the conductive solution (also called the basic solution) becomes unstable for 
a critical vertical temperature gradient beyond a certain threshold. Typically, 
instabilities in Rayleigh-Benard convection lead to roll patterns, however spiral 
waves have been also reported, in both theoretical and experimental results. 
First experimental observations on spirals |21 01 in Rayleigh-Benard convection 
were rather surprising as they were carried out for a set of external parameter 
values where parallel rolls should be stable. In experiments, spirals have been 
reported for finite Prandtl number ^ and both single and multiarmed giant spi- 
rals and spiral chaos have been observed. In |S] transitions between spiral and 
target states are reported. First theoretical studies attempting to justify spiral 
patterns in Rayleigh Benard convection were performed just in model equations 
such as Swift-Hohenberg type equations [SI [7] . These results were confirmed by 
later simulations of the full governing equations [HI Hj . 

A more general setup for Rayleigh-Benard convection than the case of uni- 
form heating, consists of a basic dynamic flow imposed by a non-zero horizontal 
temperature gradient which may be either constant or not. Numerical results 
obtained at the former setup in an annular domain address the importance of 
both vertical and horizontal temperature gradients in the development of in- 
stabilities ^1 1201 122] ■ In this work we consider a cylidrical domain with a 
localised Gaussian-like heating at the bottom. Therefore, apart from vertical 
and horizontal temperature gradients, a heat parameter related to the shape of 
the heating profile at the bottom is introduced. Our results show that many 
different types of instabilities may appear depending on the basic flow and exter- 
nal parameters. We study how instabilities depend on dimensionless parameters 
and we recover many results as those described for the classical uniform heating 
case. Giant single armed spiral waves appear and their instability threshold is 
close to those patterns refered as targets. Localised and extended spirals are 
identified and both of them are found for finite and infinite Prandtl number. 
Moreover, stationary patterns appear both with local and extended structure. 

The paper is organized as follows. Section 2 describes the physical setup, 
the general mathematical formulation of the problem in a dimensionless form, 
for basic solutions and their linear stability analysis. Section 3 explains the 
numerical method in detail and presents a test of convergence. In section 4 the 
numerical results on both, the basic and growing perturbations, are discussed 
and conclusions are presented. 

2 Formulation of the problem 

The physical set up (see figure 1) consists of a horizontal fluid layer in a cylindri- 
cal container of radius I (r coordinate) and depth d (z coordinate). The upper 
surface is open to the atmosphere and the bottom plate is rigid. At z = the 
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imposed temperature is a Gaussian profile which takes the value T^ax at r = 
and the value Tmin at the outer part {r = I). The environmental temperature 
is To. We define AT, = T^.. - To, AT,, = T„,ax - T^in and S ^ ATh/AT,. 

In the equations governing the system Ur, and Uz are the components 
of the velocity field u, T is the temperature, p is the pressure, r is the radial 
coordinate and t is the time. The magnitudes are expressed in dimensionless 
form after rescaling in the following form: r' = r/d, t' = Kt/(P , u' = du/n, p' = 
<Pp/ [poKiy] , Q = {T — To) / AT. Here r is the position vector, k the thermal 
diffusivity, h' the kinematic viscosity of the liquid and po is the mean density at 
the environment temperature Tg. After reescaling the domain f2i = [0, 1] x [0, d] 
is transformed into O2 = [0,r] x [0, 1] where T = l/d is the aspect ratio 

The system evolves according to the momentum and the mass balance equa- 
tions and to the energy conservation principle, which in dimensionless form are 
(the primes in the corresponding fields have been dropped), 

V • u = 0, (1) 

9*9 + u • ve = V^e, (2) 

dtu + (u • V) u = Pr (-Vp + V\ + ROe^) , (3) 

where the operators and fields are expressed in cylindrical coordinates and the 
Oberbeck-Bousinesq approximation has been used. Here is the unit vector 
in the z direction. The following dimensionless numbers have been introduced: 
the Prandtl number Pr= v/ k and the Rayleigh number R= gaATd'^ / kv, which 
represents the buoyant effect. In these definitions a is the thermal expansion 
coefficient and g is the gravity constant. 

We discuss now the boundary conditions (be). The top surface is flat, which 
implies the following condition on the velocity, 

= 0, on z = 1. (4) 

and free slip, 

dzUr = 0, dzU^ = 0, on z = 1. (5) 
Lateral and bottom walls arc rigid, so 

Ur ~ U0 = Mz = 0, on z = 0, (6) 
Ur = U0 = Mz = 0, on r = F. (7) 

For temperature we consider the dimensionless form of Newton's law for heat 
exchange at the surface, 

dzQ = -BQ, on z = 1, (8) 
where B is the Biot number. At the bottom a gaussian profile is imposed, 

9 = 1- 5(e^^)' - e(^"'r)'^)')/(e(^)' - 1) on z = 0, (9) 
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where /3 is a measure of the sharpness of the profile. In figure 2 several plots of 
this profile for different values of the parameters can be seen. The lateral wall 
is insulating, 

drO = on r = T. (10) 

The use of cylindrical coordinates, which are singular at r = 0, requires regu- 
larity conditions on velocity, pressure and temperature fields. In general, these 
conditions are expressed as follows [30] 

= "4>'^^ = ^<)'P on r = 0, (11) 

where and are the unit vectors in the r and cj) directions respectively. 
To summarize, the dimensionless equations contains these external parameters 
(R, r, Pr, i5, S, /3). Our aim is to describe the type of bifurcations that appear 
when some of those dimensionless quantities are varied. 



2.1 Basic state and linear stability analysis 

The horizontal temperature gradient at the bottom (i.e, 5^0) settles a sta- 
tionary convective motion which is called basic state. It is a time independent 
solution to the stationary problem obtained from equations (|1I3|I . The basic 
state is axisymmetric therefore it depends only on r — z coordinates (i.e. all 
(p derivatives arc zero). The velocity field of the basic flow is restricted to 
u = {ur,u^ — 0,Uz). Regularity conditions now become 

Ur ~ drUz = drP = 9^0 = On r = 0. (12) 

We have solved numerically the simplified equations for the basic state together 
with its boundary conditions. We use a Chebyshev collocation method with 
details given in section 3. 

The stability of the basic state is studied by perturbating it with a vector 
field depending on the r, cj) and z coordinates, in a fully 3D analysis: 



Ur{r,4>,z) = u'^{r, z) + Ur{r, z) exp{ik(l) + Xt), (13) 

u^{r,4>,z) = u''^{r,z) + Uii,{r,z)c-Kp{ik4> + Xt), (14) 

Uz{r^4>,z) = u\{r, z) + Uzir, z) cxpiikip + \t), (15) 

e(r,0,z) = Q\r,z) + Q{r,z)cyip{ik(j) + \t), (16) 

p{r,4),z) = p^{r, z) +p{r, z) exp{ik(j) + Xt). (17) 



Here the superscript b indicates the corresponding quantity in the basic state and 
the bar refers to the perturbation. We have considered Fourier mode expansions 
in the angular direction, because along it boundary conditions are periodic. 
Expressions I|13I17|) are replaced into basic equations (|1I3|I and the resulting 
system is linearized. Boundary conditions for perturbations {ur,u^,Uz,&,p) 
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are found by substituting (|13I17|) into (|4I11|I . Regularity conditions depend 
now on the wavenumbcr k: 

du, dp 
Ur ^ u^^ — — = — = — = 0, for fc = 0, (18) 
or or or 

Ur + iUtp =^ ~ & ~ p ~ 0, for k ~ I, (19) 
Mr = U0 = = e = p = 0, for k^O,l. (20) 

The resulting problem is an eigenvalue problem in A. If i?e(A) < for all 
eigenvalues the basic state is stable while if there exists a value of A such as 
i?e(A) > the basic state becomes unstable. The condition Re{X) ~ may be 
satisfied for certain values of the external parameters, (R, F, Pr, S, B, (i), which 
define the critical threshold. At the critical threshold, a stationary bifurcation 
takes place if Im{X) = while it is a Hopf bifurcation if Im{\) ^ 0. 



3 Numerical method 

We have solved both the basic state and the linear stability problem as stated 
in H1I3() . i.e. in primitive variables formulation by expanding the fields with 
Chebyshev polynomials (see Ref. |5S]). Using this technique, the problem of 
the spurious modes for pressure arises |28l I29| . which we have solved using 
the method proposed in [211, taking additional boundary conditions. They are 
obtained by the continuity equation at z = and the normal component of the 
momentum equations on r = F and z = 1, 

V • u = 0, on z 0, (21) 

p / dUr dUr dUr dUr '^\\ 

\ dt ^ dr r 86 ^ dz r ] 



^TT + ^ ^ ^ onr = F, (22) 

or Ocp 



Pr-i 



dt ^ dr r 86 ^ dz 



dv 

= -^ + Au^+Re, onz=l, (23) 



dz 

where A = r~^d / dr{rd / dr) + r^'^d'^ / dcjP' + d'^ /dz 



3.1 Basic state 

Wc have solved numerically the stationary axisymmetric version of equations 
H1I3|) together with the boundary conditions, by treating the nonlincarity with 
a Newton-like iterative method. In a first step the nonlincarity was neglected 
and a solution was found by solving the linear system: u^, u°, w^, p^ , . This 
solution was corrected by perturbation fields: u], = + u^, u\ = u'^ + u^, 
pi = p° + p and 8^ = + 0. These expressions arc introduced into the 
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equations which are Hnearized around the approach at step 0. This procedure 
is repeated with the new solution and so on. The criterion of convergence 
considered to stop the iterative procedure is that the P norm of the computed 
perturbation should be less than 10~^. 

At each step the resulting linear system for perturbations is solved by ex- 
panding any unknown perturbation field x in Chebyshev polynomials, 

L-l Af-1 

^ = E E ^TnTi{r)T„{z). (24) 

1=0 n=0 

For computational convenience the domain Q2 = [0, F] x [0, 1] is transformed 
into = [—1, 1] X [—1, 1]. This change of coordinates introduces scaling factors 
in equations and boundary conditions which are not explicitly given here. There 
are 4 x L x unknows which are determined by a collocation method. In partic- 
ular, expansions (|24fl arc replaced into the linearized, stationary, axisymmctric 
equations and boundary conditions, and those are posed at the Gauss-Lobatto 
collocation points (r^, z^), 

J-1 



L - 1 
i - 1 
N - 1 



-IjTrj , V j = O...L. (25) 
- 1 ) TT ) , V i = 0...N. (26) 



Evaluation rules arc as follows: the conveniently simplified Eqs. Ij 113(1 are eval- 
uated at nodes i = 2, ...,7V — 1, j = 2, ...,L — 1; the boundary conditions at 
z = ^1, © and lO at i = 1, j = 2, L — 1; the boundary conditions at 
z = 1, (|4I5|I and ^ a.t i ~ N, j ~ 2, ...,L — 1 ; the boundaries at r = — 1, 
ifT^ at i = 1,...,A, j ~ 1; finally the boundaries at r = 1, |7J) and ((Tn|l at 
i = 1,...,7V, j = L. The system of equations is completed with additional 
boundary conditions that eliminate spurious modes for pressure. At z = — 1 
equation (^1)) is evaluated at nodes i = 1, j — 2, L ~ 1; a,t z — 1 equation 
H23|l is evaluated at nodes i = N, j = 2, at r = 1 equation 1(22(1 is evalu- 

ated at i = 1, A^ — 1, j ~ L. With these rules the matrix associated to the 
linear algebraic system is singular, due to the fact that pressure is defined up to 
an additive constant. To fix this constant the boundary condition ((22(1 in node 
i — N — 2,j — Lis replaced by a Dirichlet condition for pressure (i.e, p = at 
i = N — 2, j = L). In this way at each iteration step a linear system of the form 
AX = B is obtained, where X is a vector containing Ax L x N unknowns and 
A is full rank matrix of order P x P with P = A x L x N . This can be easily 
solved with standard routines. 



3.2 Linear stability analysis 

The eigenvalue problem described in section 3.1 is discretized by expanding per- 
tubations ((13I17() in a truncated series of orthonormal Chebyshev polynomials 
as we did with the basic state. These expressions are replaced into the linearized 
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version of equations and boundary and regularity conditions H4I11|I . We 

use as in the previous subsection a collocation method where equations are eval- 
uated at the Gauss-Lobatto points. Evaluation rules are analogous except that 
now regularity conditions are replaced by those expressed in H18I20|I . We notice 
that for fc = 1 there are only four boundary conditions at r — —1, and there- 
fore in order to fit the unknowns with the equations we have diminished the 
expansion range of the field by one order, 

L-2 Af-1 
"■^ = X! X! a/'nTi(0Tn(2). 

;=0 ri=0 

In this way for fc = 1 we have P^4xNxL + Nx (L — 1) unknowns and 
equations while these are P = 5xixA^forfc>l and P = 4 x L x iV for fc = 0. 
For the eigenvalue problem these rules are explicitly detailed in Ref. |3D] with 
the sole difference that here equations are linearized around a non trivial basic 
state instead of around a conductive solution. 

The eigenvalue problem is then transformed into its discrete form 

Aw = \Bw (27) 

where w is a vector which contains P unknowns and A and B are PxP matrices. 
The discrete eigenvalue problem (j27|l has a finite number of eigenvalues Xi . The 
stability condition explained in previous subsection must be required now upon 
Amax where Amax = maxi?e(Ai). In Ref. |33j many details on how to solve this 
problem efficiently are given. 

As in 1201, we have carried out a test on the convergence of the method 
that let us assure the correctness of the results. Table I shows some results on 
convergence rates. We find that expansions of order 33 x 9 arc enough to ensure 
accuracy within 1%. 

4 Numerical Results 

In this section we describe numerical results found for different values of ex- 
ternal parameters. The problem formulated through Eqs. depends on 
the dimensionless numbers (R, F, Pr, S, B, (3). In particular we are interested in 
analyzing how those parameters related to heat conditions impose basic flows 
leading to a great variety of instabilities. Some of them such as giant spirals 
or targets have been already described in experiments for the case of Rayleigh- 
Benard convection with uniform heating [3 E |S] . 

Issues related to horizontal and vertical temperature differences have already 
been addressed in |2UII21ll^ but here we study for the first time the influence of 
the shape parameter /3. We restringe ourselves to large and medium P values, i.e. 
/? > 0.5. In this region the numerical method has good convergence and among 
other instabilities a giant spiral is observed. Smaller (3 values will be adressed 
in future works. We also analyze the influence of S meassuring the quotient 
between lateral and vertical temperature gradients. Figure 2a) illustrates the 
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dependence of the temperature boundary condition on (3 while figure 2b) shows 
its dependence on 5. 

4.1 Basic states 

We have solved numerically the stationary axisymmctric version of equations 

together with the boundary conditions, as explained in section 3.1. Different 
basic state solutions are obtained depending on the parameters. Figure 3a) 
displays a return flow with clearly inverted temperature gradients and a single 
roll on the velocity field. Figure 3b) shows a linear temperature field with no 
inversions on the temperature gradient and several corotative rolls. Figures 
3c) and 3d) displays intermediate states. In Figure 3c) we sec a rather vertical 
temperature gradients with slight inversions and a slight corrotativc roll. Figure 
3d) displays a return flow where the outer part of the cylinder is hotter than 
the inner one. These basic flows become unstable to different spatial structures 
as we explain next. 

4.2 Instabilities 

We have studied numerically the linear stability of the numerical basic states 
following explanations of section 3.2. We consider the external parameter R as 
control parameter. By control parameter we mean the parameter that changed 
leads to an instability while all the others are fixed. This occurs when Aniax(R.) 
changes from a negative value to a positive one as R varies. The value of Rc 
for which Amax(Rc) = is the critical value. Figure 4 displays Amax(Rc) as a 
function of the wave number k and all other parameters fixed. The eigenvalue 
with maximum real part corresponds to fc = 15 and as it is purely real, the 
bifurcation is stationary. Depending on the different parameters a great variety 
of instabilities are obtained cither stationary or oscillatory with different growing 
modes, which are analyzed next. 

4.2.1 Influence of (3 and 5 

In this section we study how the shape factor (3 and the temperature quotient 5 
both related to the bottom heat boundary condition affect to the instabilities. 
To start with we fix parameters F = 10, Pr = 0.4 and B = 0.05 and take 
R as control parameter. Figures 5a and 5b display critical values of Rc as a 
function of S for (3 = 0.5 and (3 = 5, respectively. Critical wave numbers kc and 
the corresponding growing modes are shown in these figures where void circles 
correspond to stationary bifurcations while crossed ones stand for oscillatory 
instabilities. Both figures show regions with qualitatively different behaviors 
which are separated by vertical lines. In figure 5a for 6 > 0.15 and in figure 
5b for 5 > 0.055 the critical R increases when 6 decreases. At the S limit 
values, the R threshold tends to grow so much that we have not found it. This 
suggests that there is an asymptote where the critical Rc tends to infinity. As 
S is increased, the control parameter R decreases and so does the critical wave 
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number. Comparing figures 5a and 5b, which differ in its (3 values, we notice 
more localized bifurcating structures for small values of /3. This part of the 
diagram is called region I. At different (3 the critical Rc tends to infinity at 
different 5 values. In figure 6, the curve Ci (solid line) that limits region I 
marks the position of this asymptote in the (3 ~ S plane. An example of a basic 
state in region I is presented in Figure 3a). It corresponds to a return flow. 

In figure 5b, at —0.05 < 5 < 0.055, an oscillatory bifurcation is obtained 
with wave number k = 0. This is detailed in figure 5c, where a zoom of this 
region is displayed. The critical Rc tends to infinity as 5 tends to the interval 
limits. The oscillatory mode fc = 1 is just below the instability threshold but 
very close to it, at the same parameter values. Figure 7 displays the values of 
Amax(R.c) as a function of the wave number k. The zero mode is the preferred 
one, but as the mode one is so close to it, it could be possible that both modes 
appear after the bifurcation. The mode k = corresponds to a target pattern 
(see figure 8a) and the mode A; = 1 to a giant one-armed spiral (see figure 8b). 
These plots represent contours of the growing temperature field eigenfunction 
and their structure are quite similar to those described in [3 El E El ■ The zone 
where we find these behaviours is called zone II (see figure 6). We do not find 
such instabilities for (3 = 0.5. In fact they exist only for large enough P values 
in a small interval around (5 = 0, excluding 6 = 0. This is clearly depicted 
in Figure 6 where curve C2 (dashed- dotted line) limits region II. Figure 3b) 
displays an example of a basic state in region II. It corresponds to a linear flow. 

In figure 5a for 6 < 0.15 and in figure 5b for —0.186 < 6 < —0.05, the 
bifurcations are stationary with a large wave number which increases when S 
decreases. The bifurcating structures are localized at the outer part of the 
cylinder as figures 5a and 5b display. In this range we can observe how for 
P = 0.5 the critical Rc decreases with 5 while for /3 = 5 Rc tends to grow to 
infinity at the limit 6 = —0.186. This suggests that there is an asymptote where 
the critical Rc tends to infinity. This behavior corresponds to region III in the 
f3 — 6 plane displayed in figure 6. Figure 3c) displays an example of a basic state 
in region III. 

Finally for 6 < —0.186 and /3 = 5 no instability has been found in zone IV. 
This means that in this case there exist external parameter sets at which the 
axisymmetric basic state is very robust and does not bifurcate to 3D structures. 
This regime is surrounded by curve C3 (dashed line) in the (3 — 5 diagram in 
figure 6. Figure 3d) displays an example of a basic state in zone IV. 

4.2.2 Influence of the rest of the parameters 

Aspect ratio V. Results in previous subsections are based on a rather large 
aspect ratio domain (F = 10). At smaller aspect ratii the geometry of the 
domain has more influence on instabilities. Typically for uniform heating at 
smaller aspect ratii instability thresholds are larger than those obtained for 
larger aspect ratii and critical wave numbers are smaller. Figure 9 shows how 
the instability thresholds and wave numbers change with F. It is observed that 
the standard rules just described for the case of uniform heating are mantained, 
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although now oscillatory instabilities are possible. The same type of instabilities 
as those previously described are observed. We notice that in figure 9 for F = 3 
a one-armed spiral is the most unstable mode while giant one-armed spirals 
described before, where found to be less preferred than target patterns with 
fc = 0. 

Pr number. Similar behaviors to those we have just described are also present 
for large and infinite Pr number. For instance the giant one-armed oscillatory 
spiral appears as well in a similar range of values of the parameters /3 and i5 (see 
figure 8c) and with the same conditions, the mode fc = 1 oscillatory is just below 
the instability threshold but very close to it, at the same parameter values. 

4.3 Discussion and conclusions 

We have studied how stationary and axisymmetric basic states bifurcate to dif- 
ferent 3D structures depending on the shape factor f3 and temperature ratio 
d = ATh//^Ty appearing on the temperature boundary condition. The influ- 
ence of horizontal and vertical temperature differences on instabilities has been 
addressed in numerous works [TUI [H US HHHl EOl EB EH ; in our results a 
new heat parameter is introduced related to the shape of the applied heat. We 
have found that at each /3, different 5 values define zones with qualitatively 
different instabilities. Typically at small 6 the basic state bifurcates to target 
and giant spiral wave patterns. The appearance of these structures is quite 
similar to those reported in experiments such as [3 El 01 E] for the uniform 
heating case. The reason for this similitude might be the presence of a mean 
flow in those experiments close to our basic state, that could be caused by a 
non strictly uniform heat source. We find that even for infinite Prandtl number 
spiral patterns are present. Larger negative 6 values originates very localized 
states at the outer part of the domain. Localized structures in the centre of the 
cell are found at small /? values and positive 6 values. A decreasing aspect ratio 
increases the instability threshold and decreases the critical wavenumber. This 
is in agreement with the uniform case which corresponds to (5 = 0. 

Acknowledgements 

This work was partially supported by the Research Grants MCYT (Spanish 
Government) BFM2003-02832, CCYT (JC de CastiUa-La Mancha) PAC-05- 
005-01/02, MEG (Spanish Government), MTM2004-00797, SIMUMAT S-0505- 
ESP-0158 (Gomunidad de Madrid) and by the University of Gastilla-La Man- 
cha. AMM acknowledges MGYT (Spanish Government) for a Ramon y Gajal 
Research Fellowship. 

References 

[1] H. Benard, "Les tourbillons cellulaires dans une nappe liquide. "i?ew. Gen. 
Sci. Pures Appl. 11, 1261 (1900). 



10 



[2] E. Bodenschatz, J. R. de Bruyn, G. Ahlers, and D. S. CanncU, "Transitions 
between patterns in thermal convection". Phys. Rev. Lett 67 3078, (1991). 

[3] S.W Morris, E. Bodenschatz, D. S. CanneU and G. Ahlers, "Spiral defect 
chaos in large Aspect Ratio Rayleigh-Benard convection" . Phys. Rev. Lett 
71 2026, (1993). 

[4] B. B. Plapp, D. A Egolf, E. Bodenschatz. W. Pesch, "Dynamics and se- 
lection of giant spirals in Rayleigh-Benard convection" . Phys. Rev. Lett 81 
5334 (1998). 

[5] M. Assenheimer, V. Steinberg, "Transition between spiral and target states 
in Rayleigh-Benard Convection". Nature 367 (6461) 345 (1994). 

[6] M. Bestehorn, M. Fantz, R. Fricdrich, H. Haken, C. Perez-Garcia, "Spiral 
patterns in thermal convection" . Zeitschrift fur Physik B- Condensed Matter 
88 (1): 93 (1992). 

[7] Xiao-jun Li, Hao-wen Xi, J.D. Gunton, "Dynamical properties of mul- 
tiarmcd global spirals in Rayleigh-Benard convection" , Phys. Rev. E 54 
R3105, (1996). 

[8] W. Pesch, "Complex spatiotemporal convection patterns" Chaos 6, 348 
(1996). 

[9] W. Decker, W. Pesch, and A. Weber "Spiral defect chaos in Raylcigh-Bnard 
convection" Phys. Rev. Lett 73 648, (1994). 

[10] M. K. Smith and S. H. Davis, "Instabilities of dynamic thermocapillary 
layers. 1. Convective instabihties." J. Fluid Mech. 132, 119 (1983). 

[11] A. M. Mancho, H. Herrero and J. Burguete, "Primary instabilities in con- 
vective cells due to non-uniform heating." P/iys. Rev E 56, 2916 (1997). 

[12] H. Herrero and A. M. Mancho "Influence of aspect ratio in convection due 
to non-uniform heating." Phys. Rev E 57, 7336 (1998). 

[13] R.J. Riley and G.P. Neitzel, "Instability of thermocapillary-buoyancy con- 
vection in shallow layers. Part 1. Characterization os steady and oscillatory 
instabilities." J. Fluid Mech 359 , 143 (1998). 

[14] F. Daviaud and J.M. Vince, "Traveling waves in a fluid layer subjected to 
a horizontal temperature gradient." P/i?/s. Rev. E 48, 4432 (1993). 

[15] J. Burguete, N. Mokolobwiez, F. Daviaud, N. Garnier and A. Chiffaudel, 
"Buoyant-thermocapillary instabilities in extended layers subjected to a 
horizontal temperature gradient." P/iys. Fluids 13, 2773 (2001). 

[16] J. F. Mercier and C. Normand, "Buoyant-thermocapillary instabilities of 
differentially heated liquid layers." P/i^/s. Fluids 8, 1433 (1996). 



11 



[17] M.A. Pelacho and J. Burguete, "Temperature oscillations of hydrother- 
mal waves in thermocapillary-buoyancy convection." P/ij/s. Rev. E 59, 835 
(1999). 

[18] A.M. Mancho and H. Herrero, "Instabilities in a lateraly heated liquid 
layer." Phys. Fluids 12, 1044 (2000). 

[19] A.B. Ezcrsky, A. Garcimartm, J. Burguete, H.L. Mancini and C. Pcrcz- 
Garcia, "Hydrothcrmal waves in Marangoni convection in a cylindrical con- 
tainer." F/iys. Rev. E 47, 1126 (1993). 

[20] S. Hoyas, H. Herrero and A.M. Mancho, "Thermal convection in a cylindri- 
cal annulus heated laterally. "J. Phys. A: Math and Gen. 35, 4067 (2002). 

[21] S. Hoyas, H. Herrero and A.M. Mancho, "Bifurcation diversity of dynamic 
thcrmocapillary hquid layers." P/it/s. Rev. E, 66, 057301 (2002). 

[22] S. Hoyas, A.M. Mancho, H. Herrero, N. Garnier, A. Chiffaudel. "Benard- 
Marangoni convection in a differentially heated cylindrical cavity" Phys. 
Fluids, 17, 054104 (2005). 

[23] N. Garnier and A. Chiffaudel. "Two dimensional hydrothcrmal waves in an 
extended cyHndrical vessel." Eur. Phys. J. B 19, 87 (2001). 

[24] N. Garnier and C. Normand. "Effects of curvature on hydrothcrmal waves 
instability of radial thcrmocapillary flows." C. R. Acad. Sci., Ser. IV 2 (8) 
1227 (2001). 

[25] E. Favre, L. Blumcnfeld and F. Daviaud, "Instabilities of a liquid layer 
locally heated on its free surface." P/ij/s. Fluids 9, 1473 (1997). 

[26] B.-C. Sim, A. Zebib. "Effect of free surface heat loss and rotation on tran- 
sition to oscillatory thcrmocapillary convection." Phys. Fluids 14 (1), 225 

(2002) . 

[27] B.-C. Sim, A. Zebib, D. Schwabc. "Oscillatory thcrmocapillary convection 
in open cylindrical annuli. Part 2. Simulations." J. Fluid Mech. 491, 259 

(2003) . 

[28] C. Bernardi and Y. Maday. Approximations spectrales de problemes aux 
limites elliptiques. Springer- Verlag, Berlin, (1991). 

[29] C.Canuto, MY Hussainc, A. Quartcroni, TA Zang. Spectral Methods in 
Fluid Dynamics. Springer: Berlin, 1988. 

[30] H. Herrero and A.M. Mancho, "On pressure boundary conditions for ther- 
moconvective problems." /ni. J. Numer. Meth. Fluids 39, 391 (2002). 

[31] S. Hoyas. Estudio tcorico y numcrico de un problcma dc convcccion. Thesis, 
Universidad Complutense de Madrid, 2003. 



12 



[32] N. Garnier. Ondes non-lineaires a une et deux dimensions dans une mince 
couche de fluidc. Thesis, University Paris 7 Denis Diderot, 2000. 

[33] M.C. Navarro, A. Wathcn, A. M. Mancho, H. Hcrrero. Journal of Compu- 
tational Physics submitted. 



13 



Figure captions 

Figure 1 
Physical set-up. 
Figure 2 

Temperature boundary conditions at the bottom for F ~ 10; a) for (5=1 and 
different values of /? ; b) for = 0.5 and different values of S. 
Figure 3 

a) Isotherms and velocity field of the basic state corresponding to values of the 
parameters R = 356880, F = 10, Pr = 0.4, B = 0.05, S = 0.25 and /3 = 0.5; 

b) isotherms and velocity field of the basic state corresponding to values of the 
parameters R = 19064, F = 10, Pr = 0.4, B = 0.05, 6 = 0.05 and /3 = 5; c) 
isotherms and velocity field of the basic state corresponding to values of the 
parameters R = 30377, F = 10, Pr = 0.4, B = 0.05, 6 = 0.15 and f3 = 0.5; 
d) isotherms and velocity field of the basic state corresponding to values of the 
parameters R = 10250, F = 10, Pr = 0.4, B = 0.05, S = -0.6 and (3 = 5. 
Figure 4 

Maximum real part of the growth rate A as a function of k for basic state 
at threshold Rc = 15095, the rest of the parameters arc F = 10, Pr = 0.4, 
B = 0.05, S = 0.05 and (3 = 0.5. The maximum determines the critical kc = 15. 
Void circles correspond to real eigenvalues while crossed ones stand for complex 
eigenvalues. 
Figure 5 

a) Critical R, wave number k values and growing modes as a function of S for 
P = 0.5. Remaining external parameters are F = 10, Pr = 0.4 and B = 0.05; 

b) the same for /3 = 5. Void circles correspond to stationary instabilities while 
crossed ones stand for oscillatory instabilities; c) A zoom of zones II and III of 
b). Dashed vertical lines correspond to the asymptotes. 

Figure 6 

Different regions of bifurcations in the (3 — S plane. The rest of the parameters 
are F = 10, Pr = 0.4 and B = 0.05. 
Figure 7 

Maximum real part of the growth rate A as a function of k for basic state at 
threshold Rc = 19064 corresponding to the target (fee = 0) and the giant one- 
armed spiral (fc = 1) just below the instability threshold but very close to it. 
The rest of the parameters are F = 10, Pr = 0.04, B = 0.05, 5 = 0.05 and 
(3 = 5. Void circles correspond to real eigenvalues while crossed ones stand for 
complex eigenvalues. 
Figure 8 

a) Target at the instability threshold R = 19064. The rest of the parameters 
are F = 10, Pr = 0.4, B = 0.05, 5 = 0.05, and (3 = 5]h) giant one-armed spiral 
just below the instability threshold but very close to it, at the same parameter 
values; c) giant one-armed spiral for Pr = oo, the rest of the parameters are 
R = 4674, F = 10, B = 0.2, d = 0.1 and (3 = 5. 
Figure 9 

Critical R, wave number k values and growing modes as a function of F for 
P = 0.5, Pr = 0.4, B = 0.05 and S = I. 
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Table captions 

Table I 

Critical Raylcigh numbers Rc for different orders of expansions in Chebyshev 
polynomials. The critical wave number is A:c = and the parameters are T — 10, 
Pr = 0.4, S ^ 0.05, B = 0.05 and /? = 5. 





N = 9 


iV = 11 


TV = 13 


iV = 15 


N ^17 


L 


= 19 


19069 


19070 


19070 


19070 


19070 


L 


= 25 


19055 


19056 


19056 


19056 


19056 


L 


= 29 


19065 


19066 


19066 


19066 


19066 


L 


= 33 


19063 


19064 


19064 


19064 


19064 


L 


= 37 


19063 


19064 


19064 


19064 


19064 



15 



Figure 1 



d 




16 



Figure 2 




17 



Figure 3 




18 



Figure 4 




19 



Figure 5 




20 



21 



Figure 6 




* * 



-*- 



* 



* 



— — ^s- 



/ 



IV 



— * — 




-*- ■ 




- * - 






5 = 



-0.8^ 
0.5 



1 1.5 



2.5 3 



3.5 



4.5 



22 



Figure 7 




23 




25 



